% Digital Image Correlation - Visualization of Data
%
% Author: David Gayman
% Date: 12-31-2011
%
% Instructions
%   1. Run ./correlate with option flags, as described in the
%        documentation, CORRELATE_Manual.pdf
%   2. In the internally-specified output directory, the
%        'positionField.txt' and 'strainField.txt' files will be updated
%   3. Run this script from the output directory
%   4. Visual output will be generated
%
% Interpretation of output
% %   3D Plot:
% %     -Each dot has been identified by the c++ code as a "dot"
% %     -A line connects each dot with its corresponding dot from the
% %     previous image
% %     -If a dot's location has been stored as (-1, -1), it is considered
% %     invalid and is not plotted; thus, points which have no line
% %     connecting to them signify one or more unidentified dots in previous
% %     images
% %     -Black dots are in the first image, called the reference image
% %     -Green dots have not migrated more than migration_dist_thresh from
% %     the previous image (or identified location)
% %     -Red dots have migrated excessively


%% Initialize
clc
positionField = importdata('positionField.txt');
strainField = importdata('strainField.txt');
% max_num_pts = 1000;
migration_dist_thresh = 60;

pf_ref = zeros(1,2);
pf = zeros(1,2);
num_images = (size(data,1)/max_num_pts) - 1;

neighbors = data(1:max_num_pts, :);
data = data(max_num_pts+1:end, 1:2);

neighbor_count = size(neighbors,2);

% Determine number of points
numPts = 0;
while 1 == 1
    
end


%% Plot position field for each image

% Initialize
pf_ref = zeros(1,2);
pf = zeros(1,2);
pf_old = zeros(1,2);
figure(2); clf; hold on; axis([1 850 1 num_images 1 850]); view(0,0)
set(gca, 'ZDir', 'rev'); grid on
xlabel('Column (px)'); ylabel('Image Number'); zlabel('Row (px)')

% Loop through each image i and intelligently plot points
for i = 1:num_images
    if i == 1
        % For the first image, store the actual number of points and the
        % reference position field
        for p = 1:max_num_pts
            if data(p,1) ~= -1
                num_pts = p;
            end
        end
        
        % Store reference position field
        pf_ref = data(1:num_pts, :);
        
        % Plot reference position field
        for j = 1:num_pts
            if pf_ref(j,1) ~= -1
                plot3(pf_ref(j,2), i.*ones(num_pts,1), pf_ref(j,1), 'k.')
            end
        end
        
        % Store old position field
        pf_old = pf_ref;
    else
        % For images other than reference image, update position field
        pf = data( (i-1)*max_num_pts+1:(i-1)*max_num_pts+num_pts, : );
        
        % Plot position field
        for j = 1:num_pts
            dist = sqrt( (pf(j,1) - pf_old(j,1))^2 + (pf(j,2) - pf_old(j,2))^2 );
            
            % Plot, using color to identify points that have migrated
            % excessively
            if dist > migration_dist_thresh && pf(j,1) ~= -1
                % If previous point was also valid, plot connecting line
                % and point
                if pf_old(j,1) == -1
                    
                    % Handle "dropped" points
                    for ii = i-1:-1:1
                        % Look for a valid point
                        if data((ii-1)*max_num_pts + j, 1) ~= -1
                            dist = sqrt( ( pf(j,1) - data((ii-1)*max_num_pts + j, 1) )^2 + ( pf(j,2) - data((ii-1)*max_num_pts + j, 2) )^2 );
                            % Color point correctly
                            if dist >= migration_dist_thresh
                                plot3(pf(j,2), i, pf(j,1), 'r.')
                            else
                                plot3(pf(j,2), i, pf(j,1), 'g.')
                            end
                        end
                    end
                    
                else
                    plot3([pf_old(j,2), pf(j,2)], [i-1, i], [pf_old(j,1), pf(j,1)], 'r-')
                    plot3(pf(j,2), i, pf(j,1), 'r.')
                end
            elseif pf(j,1) ~= -1
                % If previous point was also valid, plot connecting line
                if pf_old(j,1) ~= -1
                    plot3([pf_old(j,2), pf(j,2)], [i-1, i], [pf_old(j,1), pf(j,1)], 'g-')
                end
                
                % Plot point
                plot3(pf(j,2), i, pf(j,1), 'g.')
            end
        end
        
        % Store old position field
        pf_old = pf;
    end
end














%{

% Determine reference position field
p = 1;
for i = 1:size(data,1)
    if data(i,1) == -1
        break
    end
    pf_ref(p,:) = data(i,:);
    p = p + 1;
end

% Initialize figure
figure(2); clf; hold on; axis([0 850 -850 0 1 size(data,1)/max_num_pts])

% Initialize loop variables
num_pts = size(pf_ref,1);
migration_dist_thresh = 60;
dist_max = -1;
dist_min = 10000;
dist_avg = 0;
count = 0;
pf = zeros(size(pf_ref));
pf_old = pf_ref;
rt = 0;  % "reset"

% Plot position field at each time step
for i = num_pts+1:size(data,1)
    
    % Overhead logic - eliminate invalid points in the data, and update
    %   previous position field
    p = mod(i,max_num_pts);
    if p > num_pts || p <= 0
        if rt == 1
            pf_old = pf;
            pf = zeros(size(pf));
            rt = 0;
        end
        continue;
    end
    rt = 1;
    
    % Assign point to position field
    pf(p,:) = data(i,:);
    
    % If a point in the current position field is farther than
    %   migration_dist_thresh away from the point in the previous field,
    %   then display it as a different color
    dist = sqrt( (pf(p,1) - pf_old(p,1))^2 + (pf(p,2) - pf_old(p,2))^2 );
    if(dist > migration_dist_thresh && data(i,1) ~= -1)
        plot(data(i,2), -data(i,1), 'r.')
    elseif(data(i,1) ~= -1)
        plot(data(i,2), -data(i,1), 'g.')
    end
    
    % Compute average delta
    if(data(i,1) ~= -1)
        dist_max = max([dist, dist_max]);
        dist_min = min([dist, dist_min]);
        dist_avg = dist_avg + dist;
        count = count + 1;
    end
end

% Plot reference points over any inferior points
plot(pf_ref(:,2), -pf_ref(:,1), 'k.')

% Command line output
format long;
dist_avg = dist_avg/count; disp(['Migration distance'])
disp(['Average (pixels):  ', num2str(dist_avg)])
disp(['(Min, max):        ', num2str([dist_min, dist_max])])


%}

%% Plot nearest neighbors for verification
figure(3); clf; hold on; axis([1 850 1 850])
set(gca, 'YDir', 'rev'); grid on
xlabel('Column (px)'); ylabel('Image Number'); zlabel('Row (px)')
i = 1;

for p = 1:max_num_pts
    
    % Plot all points
    for pi = 1:max_num_pts
        if data((i-1)*max_num_pts + pi, 1) > 0
            plot(data((i-1)*max_num_pts + pi, 2), data((i-1)*max_num_pts + pi, 1), 'k.')
        end
    end
    
    % Plot connecting lines
    found = 1;
    for n = neighbor_count:-1:1
        n_index = neighbors(p, n) + 1;  % Value of -1 is invalid; incremented, 0 is invalid instead
        if n_index <= 0
            found = 0;
            continue;
        end
        X = [data((i-1)*max_num_pts + p, 1); data((i-1)*max_num_pts + n_index, 1)];
        Y = [data((i-1)*max_num_pts + p, 2); data((i-1)*max_num_pts + n_index, 2)];
        if (n == 1 || n == 2) && X(1) ~= -1 && X(2) ~= -1  % Plot red lines between points used for the strain calculation (currently, the first two points are used)
            plot(Y, X, 'r-', 'LineWidth', 2);
        elseif X(1) ~= -1 && X(2) ~= -1
            plot(Y, X, 'b-');
        end
    end
    
    % Pause and clear figure
    if found == 1
%        pause
        clf; hold on; axis([1 850 1 850]); grid on
    else
        break;
    end
    
end


%% Plot strain field for each time step
image_sequence = 2;
enable_contour_plot = 1;
enable_surface_plot = 0;

for i = 2:num_images
    
    disp(i)
    
    epsxx = strainData( (i-1)*max_num_pts + 1:(i-1)*max_num_pts + num_pts, 1);
    epsxy = strainData( (i-1)*max_num_pts + 1:(i-1)*max_num_pts + num_pts, 2);
    epsyx = strainData( (i-1)*max_num_pts + 1:(i-1)*max_num_pts + num_pts, 3);
    epsyy = strainData( (i-1)*max_num_pts + 1:(i-1)*max_num_pts + num_pts, 4);
    
    pf = data( (i-1)*max_num_pts+1:(i-1)*max_num_pts+num_pts, : );
    
    switch image_sequence
        case 0
            X = [pf(2,1) pf(3,1) pf(1,1) pf(6,1);
                 pf(4,1) pf(5,1) pf(7,1) pf(8,1)];
            Y = [pf(2,2) pf(3,2) pf(1,2) pf(6,2);
                 pf(4,2) pf(5,2) pf(7,2) pf(8,2)];
            Z_xx = [epsxx(2) epsxx(3) epsxx(1) epsxx(6);
                    epsxx(4) epsxx(5) epsxx(7) epsxx(8)];
            Z_xy = [epsxy(2) epsxy(3) epsxy(1) epsxy(6);
                    epsxy(4) epsxy(5) epsxy(7) epsxy(8)];
            Z_yx = [epsyx(2) epsyx(3) epsyx(1) epsyx(6);
                    epsyx(4) epsyx(5) epsyx(7) epsyx(8)];
            Z_yy = [epsyy(2) epsyy(3) epsyy(1) epsyy(6);
                    epsyy(4) epsyy(5) epsyy(7) epsyy(8)];
        case 2
            X = [pf(1,1) pf(4,1) pf(2,1) pf(3,1) pf(5,1);
                 pf(7,1) pf(10,1) pf(11,1) pf(8,1) pf(9,1);
                 pf(14,1) pf(15,1) pf(16,1) pf(13,1) pf(17,1);
                 pf(20,1) pf(23,1) pf(22,1) pf(21,1) pf(24,1);
                 pf(27,1) pf(28,1) pf(25,1) pf(29,1) pf(26,1);
                 pf(31,1) pf(34,1) pf(32,1) pf(33,1) pf(35,1)];
            Y = [pf(1,2) pf(4,2) pf(2,2) pf(3,2) pf(5,2);
                 pf(7,2) pf(10,2) pf(11,2) pf(8,2) pf(9,2);
                 pf(14,2) pf(15,2) pf(16,2) pf(13,2) pf(17,2);
                 pf(20,2) pf(23,2) pf(22,2) pf(21,2) pf(24,2);
                 pf(27,2) pf(28,2) pf(25,2) pf(29,2) pf(26,2);
                 pf(31,2) pf(34,2) pf(32,2) pf(33,2) pf(35,2)];
            Z_xx = [epsxx(1) epsxx(4) epsxx(2) epsxx(3) epsxx(5);
                    epsxx(7) epsxx(10) epsxx(11) epsxx(8) epsxx(9);
                    epsxx(14) epsxx(15) epsxx(16) epsxx(13) epsxx(17);
                    epsxx(20) epsxx(23) epsxx(22) epsxx(21) epsxx(24);
                    epsxx(27) epsxx(28) epsxx(25) epsxx(29) epsxx(26);
                    epsxx(31) epsxx(34) epsxx(32) epsxx(33) epsxx(35)];
            Z_xy = [epsxy(1) epsxy(4) epsxy(2) epsxy(3) epsxy(5);
                    epsxy(7) epsxy(10) epsxy(11) epsxy(8) epsxy(9);
                    epsxy(14) epsxy(15) epsxy(16) epsxy(13) epsxy(17);
                    epsxy(20) epsxy(23) epsxy(22) epsxy(21) epsxy(24);
                    epsxy(27) epsxy(28) epsxy(25) epsxy(29) epsxy(26);
                    epsxy(31) epsxy(34) epsxy(32) epsxy(33) epsxy(35)];
            Z_yx = [epsyx(1) epsyx(4) epsyx(2) epsyx(3) epsyx(5);
                    epsyx(7) epsyx(10) epsyx(11) epsyx(8) epsyx(9);
                    epsyx(14) epsyx(15) epsyx(16) epsyx(13) epsyx(17);
                    epsyx(20) epsyx(23) epsyx(22) epsyx(21) epsyx(24);
                    epsyx(27) epsyx(28) epsyx(25) epsyx(29) epsyx(26);
                    epsyx(31) epsyx(34) epsyx(32) epsyx(33) epsyx(35)];
            Z_yy = [epsyy(1) epsyy(4) epsyy(2) epsyy(3) epsyy(5);
                    epsyy(7) epsyy(10) epsyy(11) epsyy(8) epsyy(9);
                    epsyy(14) epsyy(15) epsyy(16) epsyy(13) epsyy(17);
                    epsyy(20) epsyy(23) epsyy(22) epsyy(21) epsyy(24);
                    epsyy(27) epsyy(28) epsyy(25) epsyy(29) epsyy(26);
                    epsyy(31) epsyy(34) epsyy(32) epsyy(33) epsyy(35)];
    end
    
    % Handle invalid strain data
    %mask = ~( (X < 0) + (Y < 0) );
    %mask = mask > 0;
    %Z_xx = Z_xx.*mask;
    
    for r = 1:size(Z_xx,1)
        for c = 1:size(Z_xx,2)
            if X(r,c) < 0 || Y(r,c) < 0
                Z_xx(r,c) = NaN;
                Z_xy(r,c) = NaN;
                Z_yx(r,c) = NaN;
                Z_yy(r,c) = NaN;
            end
        end
    end
    
    % Plot contours
    if enable_contour_plot == 1
        figure(4); clf
        subplot(2, 2, 1); hold on; axis([1 850 1 850]); set(gca, 'YDir', 'rev'); grid on; title('eps_{xx}'); xlabel('Column (px)'); ylabel('Row (px)')
        contourf(Y, X, Z_xx); colorbar
        subplot(2, 2, 2); hold on; axis([1 850 1 850]); set(gca, 'YDir', 'rev'); grid on; title('eps_{xy}'); xlabel('Column (px)'); ylabel('Row (px)')
        contourf(Y, X, Z_xy); colorbar
        subplot(2, 2, 3); hold on; axis([1 850 1 850]); set(gca, 'YDir', 'rev'); grid on; title('eps_{yx}'); xlabel('Column (px)'); ylabel('Row (px)')
        contourf(Y, X, Z_yx); colorbar
        subplot(2, 2, 4); hold on; axis([1 850 1 850]); set(gca, 'YDir', 'rev'); grid on; title('eps_{yy}'); xlabel('Column (px)'); ylabel('Row (px)')
        contourf(Y, X, Z_yy); colorbar
    end
    
    % Plot surfaces
    if enable_surface_plot == 1
        figure(5); clf
        subplot(2, 2, 1); hold on; axis([1 850 min(min(Z_xx)) max(max(Z_xx)) 1 850]); set(gca, 'ZDir', 'rev'); grid on; title('eps_{xx}')
        surf(Y, Z_xx, X); view(0, 0)
        subplot(2, 2, 2); hold on; axis([1 850 min(min(Z_xy)) max(max(Z_xy)) 1 850]); set(gca, 'ZDir', 'rev'); grid on; title('eps_{xy}')
        surf(Y, Z_xy, X); view(0, 0)
        subplot(2, 2, 3); hold on; axis([1 850 min(min(Z_yx)) max(max(Z_yx)) 1 850]); set(gca, 'ZDir', 'rev'); grid on; title('eps_{yx}')
        surf(Y, Z_yx, X); view(0, 0)
        subplot(2, 2, 4); hold on; axis([1 850 min(min(Z_yy)) max(max(Z_yy)) 1 850]); set(gca, 'ZDir', 'rev'); grid on; title('eps_{yy}')
        surf(Y, Z_yy, X); view(0, 0)
    end
    
    % Wait for user to advance
    pause
    
end
% Manually identify point order
figure(6); clf; hold on; axis([1 850 1 size(pf,1) 1 850]); view(0,0)
set(gca, 'ZDir', 'rev'); grid on
xlabel('Column (px)'); ylabel('Point Index'); zlabel('Row (px)')
plot3(pf(:,2), [1:size(pf,1)], pf(:,1), 'r.')
pause


%% Color threshold analsyis
%{
Rvalues = [16 20 55 5 59 53 40 32 45 96]./255;
Gvalues = [46 44 29 51 41 66 35 60 67 31 20 19]./255;
Bvalues = [35 20 37 64 42 62 62 43 41 59 43 14 38]./255;

figure(3); clf; plot([1:length(Rvalues)],Rvalues,'r', [1:length(Gvalues)],Gvalues,'g', [1:length(Bvalues)],Bvalues,'b')
legend('R', 'G', 'B')
disp('Mean, variance, standard deviation')
disp(['R: ', num2str([mean(Rvalues), var(Rvalues), std(Rvalues)])])
disp(['G: ', num2str([mean(Gvalues), var(Gvalues), std(Gvalues)])])
disp(['B: ', num2str([mean(Bvalues), var(Bvalues), std(Bvalues)])])
%}



%% Solve SOE
%{
x0 = sym('x0');
x1 = sym('x1');
x2 = sym('x2');
y0 = sym('y0');
y1 = sym('y1');
y2 = sym('y2');

b0 = sym('b0');
b1 = sym('b1');
b2 = sym('b2');
b3 = sym('b3');
b4 = sym('b4');
b5 = sym('b5');


A = [x0 y0 1 0 0 0;
     x1 y1 1 0 0 0;
     x2 y2 1 0 0 0;
     0 0 0 x0 y0 1;
     0 0 0 x1 y1 1;
     0 0 0 x2 y2 1];
B = [b0; b1; b2; b3; b4; b5];


disp(A\B)
%}


%% Confirms data_map works
%{
clc
validpoints=[];
for i = 1:size(neighbors,1)
    if neighbors(i,1) ~= -1
        disp(' ')
        disp([neighbors(i,1:10)])
        
        validpoints = [validpoints; neighbors(i,1:10)];
        
    end
end

for i=1:size(validpoints,1)
    for j=1:4
        disp(data( validpoints(i,j), : ))
    end
end
%}


%% Debugging large strain values
clc;

% Actual data; data1 produces good strain values, data2 really big values
data1 = [743.8630  723.5150  761.4810  804.9190  634.4790  717.2100
  738.7900  720.9240  743.6870  792.1770  641.5020  710.3600];
data2 = [743.8630  723.5150  761.4810  804.9190  634.4790  717.2100
  738.1200  709.4920  740.9430  786.1980  629.0650  716.2390];

% Calculation 1
x0 = data1(1,1);
y0 = data1(1,2);
x1 = data1(1,3);
y1 = data1(1,4);
x2 = data1(1,5);
y2 = data1(1,6);
x0_star = data1(2,1);
y0_star = data1(2,2);
x1_star = data1(2,3);
y1_star = data1(2,4);
x2_star = data1(2,5);
y2_star = data1(2,6);

b0 = x0_star - x0;
b1 = x1_star - x1;
b2 = x2_star - x2;
b3 = y0_star - y0;
b4 = y1_star - y1;
b5 = y2_star - y2;

a =  (b0*y1 - b1*y0 - b0*y2 + b2*y0 + b1*y2 - b2*y1)/(x0*y1 - x1*y0 - x0*y2 + x2*y0 + x1*y2 - x2*y1);
b = -(b0*x1 - b1*x0 - b0*x2 + b2*x0 + b1*x2 - b2*x1)/(x0*y1 - x1*y0 - x0*y2 + x2*y0 + x1*y2 - x2*y1);
c =  (b0*x1*y2 - b0*x2*y1 - b1*x0*y2 + b1*x2*y0 + b2*x0*y1 - b2*x1*y0)/(x0*y1 - x1*y0 - x0*y2 + x2*y0 + x1*y2 - x2*y1);
d =  (b3*y1 - b4*y0 - b3*y2 + b5*y0 + b4*y2 - b5*y1)/(x0*y1 - x1*y0 - x0*y2 + x2*y0 + x1*y2 - x2*y1);
e = -(b3*x1 - b4*x0 - b3*x2 + b5*x0 + b4*x2 - b5*x1)/(x0*y1 - x1*y0 - x0*y2 + x2*y0 + x1*y2 - x2*y1);
f =  (b3*x1*y2 - b3*x2*y1 - b4*x0*y2 + b4*x2*y0 + b5*x0*y1 - b5*x1*y0)/(x0*y1 - x1*y0 - x0*y2 + x2*y0 + x1*y2 - x2*y1);

%disp([a .5*(d+b); .5(b+d) e])

%% Debugging large strain values

%{
clc;


% Actual data; data1 produces good strain values, data2 really big values
data1 = [743.8630  723.5150  761.4810  804.9190  634.4790  717.2100
  738.7900  720.9240  743.6870  792.1770  641.5020  710.3600];
data2 = [743.8630  723.5150  761.4810  804.9190  634.4790  717.2100
  738.1200  709.4920  740.9430  786.1980  629.0650  716.2390];

% Calculation 1
x0 = data1(1,1);
y0 = data1(1,2);
x1 = data1(1,3);
y1 = data1(1,4);
x2 = data1(1,5);
y2 = data1(1,6);
x0_star = data1(2,1);
y0_star = data1(2,2);
x1_star = data1(2,3);
y1_star = data1(2,4);
x2_star = data1(2,5);
y2_star = data1(2,6);

b0 = x0_star - x0;
b1 = x1_star - x1;
b2 = x2_star - x2;
b3 = y0_star - y0;
b4 = y1_star - y1;
b5 = y2_star - y2;

a =  (b0*y1 - b1*y0 - b0*y2 + b2*y0 + b1*y2 - b2*y1)/(x0*y1 - x1*y0 - x0*y2 + x2*y0 + x1*y2 - x2*y1);
b = -(b0*x1 - b1*x0 - b0*x2 + b2*x0 + b1*x2 - b2*x1)/(x0*y1 - x1*y0 - x0*y2 + x2*y0 + x1*y2 - x2*y1);
c =  (b0*x1*y2 - b0*x2*y1 - b1*x0*y2 + b1*x2*y0 + b2*x0*y1 - b2*x1*y0)/(x0*y1 - x1*y0 - x0*y2 + x2*y0 + x1*y2 - x2*y1);
d =  (b3*y1 - b4*y0 - b3*y2 + b5*y0 + b4*y2 - b5*y1)/(x0*y1 - x1*y0 - x0*y2 + x2*y0 + x1*y2 - x2*y1);
e = -(b3*x1 - b4*x0 - b3*x2 + b5*x0 + b4*x2 - b5*x1)/(x0*y1 - x1*y0 - x0*y2 + x2*y0 + x1*y2 - x2*y1);
f =  (b3*x1*y2 - b3*x2*y1 - b4*x0*y2 + b4*x2*y0 + b5*x0*y1 - b5*x1*y0)/(x0*y1 - x1*y0 - x0*y2 + x2*y0 + x1*y2 - x2*y1);

disp([a .5*(d+b); .5(b+d) e])

% Calculation 2
x0 = data2(1,1);
y0 = data2(1,2);
x1 = data2(1,3);
y1 = data2(1,4);
x2 = data2(1,5);
y2 = data2(1,6);
x0_star = data2(2,1);
y0_star = data2(2,2);
x1_star = data2(2,3);
y1_star = data2(2,4);
x2_star = data2(2,5);
y2_star = data2(2,6);

b0 = x0_star - x0;
b1 = x1_star - x1;
b2 = x2_star - x2;
b3 = y0_star - y0;
b4 = y1_star - y1;
b5 = y2_star - y2;

a =  (b0*y1 - b1*y0 - b0*y2 + b2*y0 + b1*y2 - b2*y1)/(x0*y1 - x1*y0 - x0*y2 + x2*y0 + x1*y2 - x2*y1);
b = -(b0*x1 - b1*x0 - b0*x2 + b2*x0 + b1*x2 - b2*x1)/(x0*y1 - x1*y0 - x0*y2 + x2*y0 + x1*y2 - x2*y1);
c =  (b0*x1*y2 - b0*x2*y1 - b1*x0*y2 + b1*x2*y0 + b2*x0*y1 - b2*x1*y0)/(x0*y1 - x1*y0 - x0*y2 + x2*y0 + x1*y2 - x2*y1);
d =  (b3*y1 - b4*y0 - b3*y2 + b5*y0 + b4*y2 - b5*y1)/(x0*y1 - x1*y0 - x0*y2 + x2*y0 + x1*y2 - x2*y1);
e = -(b3*x1 - b4*x0 - b3*x2 + b5*x0 + b4*x2 - b5*x1)/(x0*y1 - x1*y0 - x0*y2 + x2*y0 + x1*y2 - x2*y1);
f =  (b3*x1*y2 - b3*x2*y1 - b4*x0*y2 + b4*x2*y0 + b5*x0*y1 - b5*x1*y0)/(x0*y1 - x1*y0 - x0*y2 + x2*y0 + x1*y2 - x2*y1);

disp([a .5*(d+b); .5(b+d) e])

% Calculation 2
x0 = data2(1,1);
y0 = data2(1,2);
x1 = data2(1,3);
y1 = data2(1,4);
x2 = data2(1,5);
y2 = data2(1,6);
x0_star = data2(2,1);
y0_star = data2(2,2);
x1_star = data2(2,3);
y1_star = data2(2,4);
x2_star = data2(2,5);
y2_star = data2(2,6);

b0 = x0_star - x0;
b1 = x1_star - x1;
b2 = x2_star - x2;
b3 = y0_star - y0;
b4 = y1_star - y1;
b5 = y2_star - y2;

a =  (b0*y1 - b1*y0 - b0*y2 + b2*y0 + b1*y2 - b2*y1)/(x0*y1 - x1*y0 - x0*y2 + x2*y0 + x1*y2 - x2*y1);
b = -(b0*x1 - b1*x0 - b0*x2 + b2*x0 + b1*x2 - b2*x1)/(x0*y1 - x1*y0 - x0*y2 + x2*y0 + x1*y2 - x2*y1);
c =  (b0*x1*y2 - b0*x2*y1 - b1*x0*y2 + b1*x2*y0 + b2*x0*y1 - b2*x1*y0)/(x0*y1 - x1*y0 - x0*y2 + x2*y0 + x1*y2 - x2*y1);
d =  (b3*y1 - b4*y0 - b3*y2 + b5*y0 + b4*y2 - b5*y1)/(x0*y1 - x1*y0 - x0*y2 + x2*y0 + x1*y2 - x2*y1);
e = -(b3*x1 - b4*x0 - b3*x2 + b5*x0 + b4*x2 - b5*x1)/(x0*y1 - x1*y0 - x0*y2 + x2*y0 + x1*y2 - x2*y1);
f =  (b3*x1*y2 - b3*x2*y1 - b4*x0*y2 + b4*x2*y0 + b5*x0*y1 - b5*x1*y0)/(x0*y1 - x1*y0 - x0*y2 + x2*y0 + x1*y2 - x2*y1);

disp([a .5*(d+b); .5(b+d) e])



%}

